#include "fermiqcd.h"

int main(int argc, char** argv) {
  define_base_matrices("FERMILAB");
  mdp.open_wormholes(argc,argv);
  int nt,nx;
  mdp_field_file_header header = get_info(argv[1]);
  assert(header.ndim==4);
  assert(header.box[2]==header.box[1]);
  assert(header.box[3]==header.box[1]);
  nt = header.box[0];
  nx = header.box[1];
  int nc=3;
  int L[]={nt,nx,nx,nx};
  mdp_lattice lattice(4,L,default_partitioning0,torus_topology,0,2,false);
  mdp_lattice lattice3d(3,L+1);
  mdp_matrix_field V(lattice3d,3,3);
  mdp_field<mdp_real> q(lattice3d,2);
  gauge_field U(lattice,nc);
  char filename[128];
  mdp_site x(lattice);
  mdp_site x3d(lattice3d);
  set_cold(U);
  int k,i,j,c,la;
  vector<mdp_real> p(4);

  forallsites(x3d)
    V(x3d)=1;
  
  for(int mu=0; mu<4; mu++) p[mu]=L[mu]/2;
  
  U.load(argv[1]);
  for(int t=0; t<L[0]; t++){
    forallsites(x3d){
      x.set(t,x3d(0),x3d(1),x3d(2));
      V(x3d)=V(x3d)*U(x,0);
    }
  }
  
  forallsites(x3d){
    mdp_complex z=trace(V(x3d));
    q(x3d,0)=real(z);
    q(x3d,1)=imag(z);
  }
  char realvtkname[128];
  char imagvtkname[128];
  sprintf(realvtkname,"%s.real_polyakov_line.vtk",argv[1]);
  sprintf(imagvtkname,"%s.imag_polyakov_line.vtk",argv[1]);
  
  //char filename2[128];
  //sprintf(filename2,"%s..vtk",arg[1]);
  //float tc=topological_charge_vtk(U,filename2,0);
  //mdp << "test_charge=" << tc << endl;
  q.save_vtk(realvtkname,-1,0,0,false);
  q.save_vtk(imagvtkname,-1,1,0,false);
  
  mdp.close_wormholes();
  return 0;
}

